Matter-wave solitons with a periodic, piecewise-constant nonlinearity 
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Motivated by recent proposals of "collisionally inhomogeneous" Bose-Einstein condensates 
(BECs), which have a spatiaUy modulated scattering length, we study the existence and stability 
properties of bright and dark matter-wave solitons of a BEC characterized by a periodic, piecewise- 
' constant scattering length. We use a "stitching" approach to analytically approximate the pertinent 

solutions of the underlying nonlinear Schrodinger equation by matching the wavefunction and its 
derivatives at the interfaces of the nonlinearity coefficient. To accurately quantify the stability of 
\^ • bright and dark solitons, we adapt general tools from the theory of perturbed Hamiltonian systems. 

^SJ ' We show that solitons can only exist at the centers of the constant regions of the piecewise-constant 

nonlinearity. We find both stable and unstable configurations for bright solitons and show that all 
dark solitons are unstable, with different instability mechanisms that depend on the soliton location. 
We corroborate our analytical results with numerical computations. 



INTRODUCTION 



In the past few years, the study of solitary wave structures in atomic physics has received considerable attention, 
fNj predominantly because of impressive advances in the field of Bose-Einstein condensation [l|, . These investigations 
K*" ' have largely been motivated by the comprehensive mean-field description of atomic Bose-Einstein condensates (BECs) 
, by the Gross-Pitaevskii (CP) equation, a variant of the nonlinear Schrodinger (NLS) equation in which the nonlinearity, 
whose strength is proportional to the s-wave scattering length, is introduced by the interatomic interactions. This 
remarkable feature has allowed both theoretical investigations and experimental observations of bright [3, 0, 
dark 0, 0, H, 0, and gap matter- wave solitons in BECs (see also the recent review [IH). Such solitons have 
' been studied in detail in the literature in the presence of various external potentials, including harmonic traps and 

J periodic lattices. This has yielded tremendous insights into a large variety of interesting phenomena, including Bloch 

^ oscillations, Landau-Zener tunneling, modulational ("dynamical") instabilities, gap excitations, and more (see the 
reviews [HI, [H, [H, [13, and references therein). The types of matter- wave solitons that arise in a given situation 
depends on the nature of the interatomic interactions which, in turn, are characterized by the sign of the scattering 
length. Specifically, bright solitons arise in BECs with attractive interactions and negative scattering length, whereas 
^ ■ dark and gap solitons emerge in the case of repulsive interactions and positive scattering length. Gap solitons also 
require the presence of an optical lattice or superlattice potential. 

During the past decade of BEC research, a number of tools have been developed to control and manipulate matter 
waves. For example, they can be influenced, processed, and shaped using static (homogeneous and inhomogeneous) 
electric and magnetic fields [l^ : optical devices [l3l , and near-field radio- frequency devices [l^ . One can manipulate 
not only a BECs external (trapping) potential but also the interactions among the atoms that are responsible for 
the nonlinear properties and dynamics of the matter waves. The interaction among ultracold atoms can be adjusted 
experimentally by employing either magnetic [l^, [13] o r op tical Feshbach resonances [2l| in a very broad range. 
Additionally, so-called "confinement induced resonances" [22l.[23j allow one to vary the effective one-dimensional (ID) 
coupling constant of quasi-lD systems by adjusting the transversal confinement length. One can consequently vary 
the external potential while independently and simultaneously changing the strength of the nonlinearity by tuning 
the interatomic interactions. 

The manipulation of BECs using Feshbach resonances has propelled a significant number of investigations that 
have been subsequently refined substantially. Experimental achievements include the formation of bright matter- wave 
solitons and soliton trains for ^Li [sj, [3| and ^^Rb [^ atoms by tuning the interatomic interaction within a stable BEC 
from repulsive to attractive, the formation of molecular condensates [IJ], and the probing of the BEC-BCS crossover 
[1^. Moreover, theoretical studies have predicted that a time-dependent modulation of the scattering length can be 
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used to stabilize attractive two-dimensional BECs against collapse [1^ or create robust matter-wave breathers in ID 
BECs [27|] . We remark in passing that relevant contexts involving dependence of the nonlinearity coefficient on the 
evolution variable have also been of interest in nonlinear optics in describing pulse propagation in an optical fiber 
in the presence of periodically varying dispersion and nonlinearity (28l . [29j (see also the references therein) . In that 
context, a popular approach has been to look for an averaged (over the fast variation scale) equation with constant 
coefficients in which solitary- wave solutions could be sought. This approach has also been useful in the context of 
BECs in describing the averaged properties of solitary waves in the presence of a periodic, time-dependent modulation 
of the scattering length [30, 31, 32]. 

In addition to the aforementioned studies involving temporal variations of the interaction strength, it has been 
recently found that atomic matter waves exhibit novel features under the influence of a spatially varying scattering 
length and, consequently, a spatially varying nonlinearity. The resulting so-called "collisionally inhomogeneous" envi- 
ronment provides a variety of interesting and previously unexplored dynamical phenomena and potential applications, 
including adiabatic compression of matter waves (ssl . [33. atomic soliton emission and atom lasers [351] , enhancement 
of the transmittivity of matter waves through barriers [361 l37l | , dynamical trapping of matter- wave solitons [sH] , stable 
condensates exhibiting both attractive and repulsive interatomic interactions 38] , and the delocalization transition of 
matter waves |3&1 . Particular inhomogeneous frameworks that have been investigated include linear [ssl [s^ , parabolic 
(40| . random [ilf, periodic [sol. li^. lisf . and localized (step-like) f35l. H . l45ll sp atial variations. There have also been a 
number of detailed mathematical studies [i^, [13, E^. In particular, Refs. j46[[47j examined the effects of a "nonlinear 
lattice potential" (i.e., a spatially periodic nonlinearity coefficient) on the stability/instability of solitary waves and 
elucidated the interplay between drift and diffraction/blow- up instabilities. More recently, the interplay of nonlinear 
and linear potentials has been examined in both continuum [1^ and discrete [50| settings. 

Motivated by the above studies, as well as the suggestion of a piecewise-constant inhomogeneous scattering length 
in Refs. (35l.l45|. we investigate in the present work matter- wave solitons in the framework of the NLS equation with a 
periodic, piecewise-constant nonlinearity coefficient. Preparing such a setup experimentally is challenging but possible 
in principle. For example, it could be achieved by integrating optical and magnetic traps on an atom chip |16j . The 
underlying idea is as follows: An on-chip fiber-optic configuration generates a microscopic dipole waveguide with a 
strong transversal confinement |5ll| . This waveguide is augmented by a periodic array of either permanent magnets 
[5^ or current-carrying wires that create alternating zones of approximately zero and constant magnetic fields 
in the longitudinal direction. Tuning the magnetic field (via, for example, the currents in the wires) such that the 
resulting field strength leads to an atomic scattering process close to a Feshbach resonance allows one to achieve a 
substantial variation of the atomic scattering behavior in the form of an alternating, approximately piecewise-constant 
scattering length. Refs. [H, considered a localized region with positive scattering length (the rest of the BEC 
was noninteracting, providing a substantial linear regime). Here, however, we consider a scattering length given by 
a piecewise-constant function that switches periodically between either two positive or two negative values. The 
importance of this particular fully nonlinear setting is that it allows us to analytically investigate the existence and 
stability of matter-wave solitons in each region (where the nonlinearity coefficient is constant) and to subsequently 
"stitch" the constructed (local) solutions at the jump points of the piecewise-constant nonlinearity function. Our 
analysis is similar in spirit to the approach used in a situation with both linear and nonlinear, piecewise-constant 
periodicity in Ref. [53]. We employ techniques from the general theory of perturbed Hamiltonian systems, starting 
from the unperturbed limit of the completely integrable ID NLS equation, in order to analytically address the stability 
of both bright 54] and dark [55 ] matter-wave solitons. In order to quantify the accuracy of our theoretical predictions, 
we subsequently compare these findings with direct numerical computations. 

Our presentation is organized as follows. In Section|lTl we present our setup and the pertinent GP model. We study 
bright and dark matter-wave solitons, respectively, in Sections IIIII and IIVI Finally, in Section|Vl we briefly summarize 
our findings and present our conclusions. 



II. SETUP AND MODEL 



We consider a "cigar-shaped" condensate that is elongated along the z direction and strongly confined in the 
transverse {x and y) directions . Assuming that this condensate is collisionally inhomogeneous, characterized by a 
spatially varying scattering length a(z), we follow the standard approach 56] of averaging the three-dimensional GP 
equation in the transverse plane. We thereby arrive at the following ID GP equation for the longitudinal wavefunction 



iHdt'ip{z, t) = 



2m 



Tp{z,t) 



(1) 



where m is the atomic mass and uj± is the transverse confinenent frequency. 



3 



Measuring time and space in units of uJj^^ and transverse oscillator length lj_ ^ ^yh/mLu±, respectively, we reduce 
the GP in Eq. ^ to the dimensionless form, 



idt(f){z,t) = 



-\d',+g{z)\<l>{z,t)\' 



{z,t), (2) 



where = y^2\ao\tlj{z) and g{z) = a{z)/\ao\ is negative (positive) for attractive (repulsive) interatomic interactions. 
In the above expressions, og denotes the value of the scattering length of the corresponding coUisionally homogeneous 
system, for which g{z) = go = const. Note that Eq. ^ has the form of a ID NLS equation. It possesses two integrals 
of motion: the normalized number of atoms N = |0|^ and the Hamiltonian. The former is connected to the 
unnormalized number of atoms A/" through the equation J\f = {a±/2\ao\)N . 

We take the nonlinearity coefficient g{z) in Eq. ([2]) to be a piecewise-constant function of space: 

g(z) = |-90' ifO<^^°^(^'^;<^i . (3) 
^ I 51 , if Li < mod(z, L) < L ^ ^ 

Note that the nonlinearity coefhcient in Eq. Q may also be expressed as g{z) = go + Ag(z), where 

n— + 00 

Ag{z) = Ago ^ {Q {z - [nL + L^]) ~ Q {z - [n + 1]L)} , (4) 

n— — oo 

where Ago = gi ~ go and is the Heaviside step function. 

In our numerical investigations (see the discussion below), we adopt values of the parameters that are typical for 
experimental setups. In particular, we consider a quasi-lD trap with a transverse confining frequency ujj_ — 2Tr x 1400 
Hz, which fixes the temporal unit to 0.1 ms. We assume that the dimensionless chemical potential fi is of order 0(1). 
For an (attractive) ^Li condensate, this choice yields a spatial unit of 1/im and a number of atoms of « 1100. 
For the (repulsive) ®^Rb (respectively, ^"^Na) condensate, this yields a spatial unit of 0.3/xm (2.2/xm) and a number 
of atoms of « 1200 (16000). Moreover, we assume that the nonlinearity function g{z) is characterized by changes 
Ag of the order of 50% of the typical values of the scattering lengths, and that the normalized periodicity L is of the 
order 0(1) (i.e., it is about a few microns). This choice is consistent with typical Feshbach resonances that occur in 
Li, Rb, and Na condensates: Changing the magnetic field by a factor of two on a micron scale in the neighborhood of 
a Feshbach resonance requires magnetic field gradients in the range 1 kG/cm to 100 kG/cm, which are of moderate 
size for atom chips. The lower limit is also compatible with the gradients used in macroscopic trap setups. 

Below, we will quantify the existence and stability properties of bright and dark matter- wave solitons (for negative 
and positive nonlinearity coefficients g{z), respectively) as a function of the inhomogeneity strength Ago- 



III. BRIGHT MATTER- WAVE SOLITONS 



A. Perturbation analysis 

Let us first consider an attractive BEG characterized by a negative nonlinearity coefficient in Eq. ^ , so that go < 
and |A(7(z)| < \go\- In the special case of a homogeneous nonlinearity, scaled such that its coefRcient is go ~ ~1, 
Eq. ([2]) possesses a bright matter- wave soliton solution of the form 

(f>bs {z,t) = r] sech[?7(z - zo)] exp(i^i) , (5) 

where rj is the soliton's amplitude and inverse width, zq is the position of its center, and /i — —(1/2)77^ is its effective 
chemical potential. We then consider the spatial inhomogeneity as a perturbation of strength e = Ago- When e 0, 
the integrability of the original NLS in Eq. ^ is lost, but the system is still Hamiltonian. Its Hamiltonian functional 
is given hy H — Ho + tHi, where the unperturbed {Ho) and perturbed {Hi) parts of the Hamiltonian are given by 

-(|a.0l'-|(/)|^)dz, (6) 

In accord with the general perturbation theory for Hamiltonian systems developed in Ref. [H^l, the condition for the 
existence of the solution ^ under the aforementioned perturbation is that it remains an extremum of the perturbation 
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Hamiltonian Hi. In the present case, it is easy to see by considering the dependence of Hi on zq that this condition is 
satisfied provided is either at the center of a region with nonlinearity coefficient or (by symmetry) at the center 
of a region with nonhnearity coefficient gi. Therefore, we can identify stationary soUtary wave solutions centered at 
these points in the presence of the periodic, piecewise-constant nonhnearity. 

Let's consider the stabihty of these bright matter- wave sohtons. As is weh-known, the unperturbed NLS equation is 
invariant under both translation and phase/gauge transformations. If the perturbation is sufficiently small, instability 
can arise in the case of attractive nonlinearity only via the perturbation-induced breaking of one of these symmetries. 
(Note that the continuous spectrum is bounded away from the origin and for small e cannot give rise to an instability 
[5J|.) Furthermore, it is clear that in the presence of a spatially-dependent nonlinearity, the U(l) symmetry (i.e., 
the phase invariance) of the equation is preserved, so the stability ultimately should depend on the location of the 
translational eigenvalue. The eigenvalues A (and the corresponding eigenfrequencies cj, which satisfy uj"^ = — A^) can 
be obtained within the framework of the perturbation analysis of Hamiltonian systems in Ref . [s^l (see also Ref. [40] ) 
according to the equation 



det (eM - uj^D) = , (8) 



where the matrices M and D are given by 




/ _S_ &Hi_ d<pts 

M = ( 9^0 \ <5r ' 02=0 / ^ I , (9) 

In the above equations, the star * denotes complex conjugation, (,) denotes the inner product, and SHi/dcj)* is the 
functional (Frechet) derivative. 

The nonzero elements of the matrices M and D (i.e., mn, dn, and ^22) can then be calculated directly to obtain 

Tl — + 00 

mii=-ff ^ {sinh[r/([n + l]i-zo)]sech^[77([n+l]L-zo)] (11) 

n— — 00 

— sinh [rj [Li + nL — Zq)] sech^ [77 [Li + nL — Zq)]} i (12) 

where dn = rj and ^22 = ^^/v- Equation ([8]) allows one to calculate the translational eigenfrequencies to leading 
order in e, yielding 

emu — Lj'^rj = . (13) 

The sum in Eq. ([TT|) can then be calculated numerically. Note that this sum converges rapidly because its component 
functions decay exponentially, and only three terms (n = —1, 0, and -I- 1) are needed to obtain a value that is accurate 
up to ten digits. Below we will directly compare the analytical prediction provided by Eq. (|13p to the numerical 
results for different values of A(7o (and of the soliton's chemical potential /i) obtained by computations. 



B. "Stitching" of bright soliton solutions 



We can use specific features of this particular system to go beyond the above theory, which is valid for general 
perturbed Hamiltonian systems, and acquire a quantitative handle on the perturbed soliton profile. In particular, 
we use the fact that one can explicitly construct the solution of the steady state equation for each of the regions in 
which the nonlinear coefRcient is constant. For sufficiently small e, we expect the solutions to deviate little from the 
solution (/)bs- We then construct an analytical approximation to the perturbed solution by "stitching" together the 
functions that would be solutions to the respective homogeneous NLS equations, with the nonlinearity coefficients 
taking different constant values in different spatial regions. To do this, we require that the wavefunction cj) and its 
first spatial derivative dz4> be continuous across the boundaries between regions of different g. 

Let us consider, in particular, the bright soliton solution of Eq. ([2]) centered at Zm with amplitude rjm in a region 
with nonlinearity coefRcient gm'- 

(pmiz) = ^^sech[7y,„(z - z,„)] . (14) 

V 9m 
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Starting with a known value of, say, rj^ and Zm, we then need to calculate all other rjj and Zj {j ^ m) as a function 
of the known ones. Applying the aforementioned continuity conditions, we obtain 

sech[?7„(Z„i - 2„j)] = ^'"'^^ scch[ri,n+i{Z,n - Zm+i)] , (15) 



/9m ^/9m+l 



„2 



-^^sech[77,„(Z„ - z,„)] tanh[77„(Z„i - Zm)] = — '^"^^ sech[rirn+i{Zrn - Zm+i)] 

X tanh[?7,„+i(Z„ - z„+i)] , (16) 

where Zm is the coordinate of the interface between regions m and m + 1. Using Eq. (fT5|l . it is then possible to 
simplify Eq. (fT6|l to get 

Vm+iiZm - Zm+i) = arctauh \ ^"^ tanh [77,„(Z„ - z^)] \ , 

Vl-rn+l J 

from which we obtain z^+i once 77^+1 is known. We can then proceed iteratively to obtain (2^1+2, ?7m+2), etc. Note, 
however, that this method is not exact, as can be inferred by close inspection. In particular, Eqs. P^ - (|16p possess two 
unknowns, Zm+i and rjm+i, provided Zm and rjm are given. However, in principle, the quantity rjm+i is determined 
by the coefficient gi of the piecewise-constant nonlinearity (and the chemical potential /x). Hence, the system is, in 
fact, overdetermined. Nevertheless, as long as /x is not too small (so that the width of the soliton is smaller than 
or comparable to the size of a step in the periodic nonlinearity (see Fig. [1} , the above stitching procedure is fairly 
accurate, as will be discussed below when we present the results of our numerical simulations. 



C. Numerical Results 



We begin the description of our numerical findings by illustrating the solution obtained by numerically solving the 
standing wave problem corresponding to Eq. ([2]) via a fixed point (Newton-Raphson) iteration. In Fig. [1] we show an 
example of both a stable bright matter- wave soliton solution (top left) and an unstable one (bottom left). We also 
present (in the right panels of the figure) the corresponding spectral planes (Re(u;), Im(w)) of the eigenfrequencies 
oj = Re(w)-|-ilm(a;) that result from a Bogoliubov-de Gennes analysis of the NLS equation around the soliton solution. 
The absence of eigenfrequencies with nonzero imaginary part indicates stability, whereas the presence of imaginary 
or complex eigenfrequencies is a signature of linear instability. The stability features observed are reminiscent of the 
ones found in earlier works, such as Refs. [i^, 113, H^. Namely, the solution is stable when centered at a maximum 
of the magnitude (in absolute value) of the piecewise-constant, periodic nonlinearity, but it is linearly unstable when 
centered at a minimum. In Fig. [2 we also show the solutions that we obtained semi-analytically via the "stitching" 
procedure illustrated above. Observe the excellent agreement between those solutions and those computed using fully 
numerical simulations. In fact, the two solutions are virtually indistinguishable in this case and, more generally, in 
all situations with sufhciently large \^\. 

Figure [5] depicts the evolution of bright matter- wave solitons with initial configurations of the forms presented in 
Fig. [1] In particular, we present the stable case (which has parameter values /i — —1.0 and At^o — +0.5) in the top 
panel. Observe that the bright soliton remains unchanged during the evolution, as predicted by the linear stability 
analysis. On the other hand, we show the unstable case (with par ameter values /x = —1.0 and Aijo = —0.5) in the 
bottom panel. Observe that for t > 80, the soliton drifts [4^l47ll49j and starts oscillating around a stable fixed point 
(the center of a neighboring maximum of the magnitude of the periodic potential) as a result of the existence of the 
imaginary eigenfrequency. 

In Fig. [3l we show a conclusive diagram that summarizes our results for the linear stability of the bright-soliton 
solutions we constructed. The plot shows the relevant translational (squared) eigenfrequency as a function of the 
strength of the perturbation Aqq for different values of the frequency /i. It is clear that the corresponding squared 
eigenfrequency is essentially linear in Ago, giving rise to instability for Ago < and stability for At^o > 0. In all 
cases, we observe very good agreement between the theoretical prediction and the numerical results, especially when 
Ago is sufficiently small (for which we expect the theory to be most accurate). 
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FIG. 1: (Color online) Stitched bright-soliton solutions versus fully numerical bright-soliton solutions of the NLS equation. 
The top left panel shows the wave profile and the "stitched" solution for the same parameter values {fi — —1.0 and Ag > 0). 
The top right panel gives the corresponding spectral plane (Re(cj), Im(cj)) of the eigenfrequencies u) = Re((.<j) + ilm(tj). In the 
bottom panels, we show the same plots for a case with ^ — —1 and Ag < 0. The dashed curves in the left panels show the 
spatially dependent nonlinearity coefficients. The bright soliton depicted in the top panels is stable, whereas that in the bottom 
panels is unstable. Notice that for these parameters the stitched solutions is practically indistinguishable from the numerically 
exact one. 
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FIG. 2: (Color online) Time evolution of the bright matter-wave soliton for the two cases shown in Fig. [T] The top panel 
shows the contour plot of |(^(2:,f)| for parameter values ^ = —1.0 and Ago = +0.5, whereas the bottom one shows the case of 
/i = —1.0 and Ago = —0.5. The bright soliton depicted in the top panel remains unchanged during the evolution, whereas that 
in the bottom panel drifts after t > 80 and oscillates around a stable equilibrium point (i.e., around a maximum of |(;(2)|). 
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FIG. 3: (Color online) Numerical results (points) and theoretical predictions resulting from the perturbation analysis (curves, 
which are labeled "the" in the inset) for the translational squared eigenfrequency lo^ as a function of perturbation strength 
Ago for different values of the effective chemical potential of the bright matter- wave soliton. 



IV. DARK MATTER- WAVE SOLITONS 



A. Perturbation analysis 

We now consider repulsive BECs, which are characterized by a positive nonlinearity coefficient in Eq. ([2]). That is, 
f/o > 0, with |A5(z)| < go as in the previous case. In the special case of a homogeneous nonlinearity, scaled so that 
its coefficient is go = +1, Eq- (E)) possesses a dark (black) matter- wave soliton solution of the form 

(l)ds{z,t) = 77tanh[?7(z - zo)] exp(i/zi) , (17) 

where r] is the soliton's amplitude and inverse width, zq is its center position, and ^ = r]^ is its effective chemical 
potential. In the case under consideration, the solution (jl7p has nonvanishing boundary conditions at infinity, with 
the asymptotic limits (f>ds{z) — > ±77 as z — > ±c». Consequently, one can cast the problem of the piecewise-constant 
nonlinearity in the form of a perturbation, reminiscent of what was done for linear external potentials in Ref. [ssl ]. 
The perturbed part of the Hamiltonian can be written formally as 




We remark that for the perturbed solutions (see the relevant figures in the numerical section below), it is not obvious 
that this integral converges. However, for the unperturbed solution (for which the solvability condition under the 
effect of perturbation is applied), it is certainly convergent. 

Following Ref. [ss'l, one can show in a straightforward manner that for the solution to survive under the effect of 
perturbation (using the perturbation parameter e = Ago), the function 

M'{s) = \ ^(z) [ry4 - _ ,)] (19) 
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must vanish. As was the case for the bright matter-wave sohtons we analyzed in the previous section, this condition 
impHes that the center of the dark matter-wave soUton must be located at the center of either a region in which the 
nonlinearity coefficient is go or at the center of one in which it is gi . 

The stability of the dark-soliton solution we constructed depends on the sign of the second derivative of M{s) 
evaluated at the root sq of the first derivative p9)) . Based on the analysis of Ref. [55|], an instability will be present 
with one imaginary eigenfrequency pair if eAf"(so) < and with exactly one complex eigenfrequency quartet if 
eM"(so) > 0. As was the case for bright solitons (but with some important differences, which we discuss below), 
this instability is dictated by the translational eigenvalue, which needs to bifurcate from the origin as soon as the 
perturbation is present. For eM"{so) < 0, the relevant eigenfrequency pair moves along the imaginary axis, leading 
to an immediate instability associated with exponential growth of a perturbation along the associated eigendirection. 
For eM"(so) > 0, on the other hand, the eigenfrequency moves in principle along the real axis, which, however, is 
filled with continuous spectrum. Consequently, as a result of the opposite signature (to these eigenfrequencies) of the 
translation mode, the eigenfrequency exits as a complex quartet, signalling the presence of an oscillatory instability. 
In fact, following Ref. [5a|, we find that the relevant eigenfrequencies are determined by the characteristic equation 



(20) 



where the perturbation is given by Eq. ([4]) and we recall that the eigenvalues A are related to the eigenfrequencies uj 
through — —uP'. 

Because the two roots of M'(s) are Sqi = "-^ + -^i/2 and so2 = "-^ + Li + one can evaluate M" [sq) explicitly 
to obtain 



M"(soi) = V XI 1 



77 



tanh'' 



(21) 



and a similar expression for M"(so2)- Combining Eqs. and (PT|) yields a prediction for the location of the relevant 
translational eigenfrequency. To our knowledge, this is the first time that this has been done (albeit formally) for the 
case of a dark soliton in the presence of a spatially inhomogeneous nonlinearity. We now turn to numerical simulations 
in order to examine the accuracy of our analytical results. 



B. Numerical Results 



In Fig. m we show the profile and spectral plane of the dark matter-wave soliton when it is centered at min[(7(z)] 
(top panels) and max[(7(z)] (bottom panels). The former case is associated with the scenario in which eM"(so) > 0, 
and its corresponding spectral plane has an eigenvalue quartet (corresponding to an oscillatory instability) . The latter 
case is characterized by eAf "(sq) < 0, which results in a real eigenvalue (imaginary eigenfrequency) pair in the spectral 
plane. 

Figure [5] depicts the evolution of dark solitons with the two initial configurations shown in Fig. [4l In both cases, 
/i — —1 and |Ago| = 0.5. In the top panel, the wave is centered at min[(7(z)], and in the bottom panel it is centered 
at max[f/(2;)]. Unlike their bright counterparts, dark solitons become mobile in both cases. This mobility develops 
earlier for the soliton in the bottom panel because its eigenfrequency has a much larger imaginary part of than that of 
the soliton in the top panel. Furthermore, in the top configuration, the instability is initiated through an oscillatory 
phase that stems from its complex eigenfrequency quartet. This feature is absent in the case of a wave centered at 
max[(7(z)], as is expected when there is an imaginary eigenfrequency pair (which corresponds to a purely exponential 
instability). It is interesting to note that in both cases, as the soliton becomes mobile over the periodically nonlinear 
terrain, it radiates (a process that has been discussed elsewhere, although for a linear periodic potential; see, for 
example, Refs. and references therein) and is eventually destroyed. 

For the min[(7(z)]-centered soliton, we compare in the bottom panel of Fig. [S]the imaginary part of the eigenfrequency 
quartet (which corresponds to the real part of the eigenvalue quartet and yields the growth rate of the instability) to 
the theoretical predictions of E qs. (P(I)l - (PT|) . Although finite-size effects similar to the ones reported in Ref. ^] (see 
the detailed discussion of Ref. [60j) are present, the agreement between the theoretical prediction and the numerical 
result is good (and becomes better as the domain size increases). We show three numerical cases and observe that 
the agreement with the theoretical prediction improves (as expected) as the density of the numerical phonon band of 
eigenvalues increases (for larger domain length V). See Ref. [60] for details on this finite size effect. We also show the 
real part of the corresponding eigenfrequency for different domain lengths, for which the change in domain size does 
not play such a significant role. As shown in Fig. [71 we observe a similar level of agreement (but without the finite 
size effects, because the relevant eigenfrequency is now imaginary) for a max[5(z)]-centered dark soliton (for which 
eM"(so) < 0). 
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V. CONCLUSIONS AND FUTURE DIRECTIONS 

We have studied the existence and stabiUty of bright and dark matter-wave sohtons in Bose-Einstein condensates 
with spatiaUy-dependent, periodic scattering lengths. Our analysis was based on an analytically tractable model 
given by a Gross-Pitaevskii equation with a periodic, piecewise-constant, spatially inhomogeneous nonlinearity. In 
particular, we used techniques from the theory of perturbed Hamiltonian dynamical systems in order to obtain 
conditions for the "persistence" of the matter-wave solitons and analyze their linear stability. We thereby found that 
solitons must be centered in one of the constant regions of the piecewise-constant coefficient of the nonlinearity (which 
is directly proportional to the atomic scattering length) . We have shown that bright-soliton solutions are stable when 
localized in regions of maximal (absolute) nonlinearity and unstable when localized in regions of minimal (absolute) 
nonlinearity. This is consistent with the findings of earlier works on the NLS equation with periodic potentials (see, 
for example, Refs. [H, \^ E^). In this situation, we also presented an approximate analytical technique (called 
"stitching" ) that allowed us to match the solution in a semi-analytical fashion at the interfaces where the nonlinearity 
coefficient changes in order to obtain an accurate profile of the soliton in the presence of the spatial inhomogeneity. 
We applied a similar approach to dark matter-wave solitons (in that the translational eigenvalue is still responsible 
for the configurational stability) to show a rather different final result. Specifically, we showed that dark solitons 
centered at min[(7(z)] feature a complex eigenvalue quartet, whereas ones centered at max[g(z)], have a pair of 
unstable real eigenvalues. Therefore, both cases are unstable, but the former instability is oscillatory and the latter 
is purely exponential (typically with a larger growth rate). In all cases, we corroborated both the existence and the 
linear stability results with corresponding numerical computations, which were in good agreement with our analytical 
findings. 

It would be interesting to extend the present results to more complicated configurations. One interesting example 
along the lines of Ref . [4^ might be to examine the competition of a piecewise- nonlinear and a piecewise- linear potential 
in one dimension. Another possibility would be to investigate the stitching between spatially extended solutions given 
by elliptic functions. Perhaps a more appealing example for a generalization of the present context would be to 
examine piecewise-constant nonlinearities in higher dimensions (with, for example, a piecewise-constant radial form 
or a piecewise-constant, square-lattice form) and determine the stability of vortex-like and similar structures in the 
defocusing NLS. In the focusing case, it would also be interesting to determine whether such a nonlinearity could 
stabilize solitary waves (which would be given piecewise generalizations of the Townes soliton) against collapse. 
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FIG. 4: (Color online) Numerical dark matter- wave soliton solutions of the NLS Eq. ((2)). The top and bottom left panels 
show the wave profiles (solid curves) in the presence of a periodic, piecewise-constant coefficient in front of the nonlinearity 
(dashed curves). In each case, the chemical potential is /i = 1. The soliton is centered at the min[(7(2)] in the top panels and at 
max[(;(z)] in the bottom ones. The right panels show the eigenfrequency spectrum (similar to the right panels of Fig. [T]). The 
top panel illustrates an oscillatory growth due to a complex quartet of eigenfrequencies, and the bottom shows an exponential 
growth due to an imaginary eigenfrequency pair. 
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FIG. 5: (Color online) Time evolution of the dark soliton for the two configurations shown in the left panels of Fig. |4] The 
dark soliton shown in the top panel remains unchanged for a longer time, showing evidence of an eigenfrequency with a smaller 
imaginary part. Furthermore, the former instability is oscillatory in its initial stages of development (indicating its association 
with a complex eigenfrequency quartet), whereas the latter is not (indicating its connection to an imaginary eigenfrequency). 
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FIG. 6: (Color online) Real (top panel) and imaginary (bottom panel) parts of the eigenfrequency as a function of perturbation 
strength Ago for a BEC with a chemical potential of = 1.0 which has a soliton centered at min[(;(2)]. The points give 
the numerical results and the solid curve gives the theoretical prediction for the real and imaginary parts of the complex 
eigenfrequency quartet. In both panels, we present the results of numerical simulations for three different domain sizes L in 
order to illustrate the role of the finite size of the computational lattice, which is especially evident in the instability growth 
rate represented by Im{u)). 




FIG. 7: (Color online) Imaginary part of the eigenfrequency as a function of the perturbation strength Ago in for a BEC with 
chemical potential jj, — 1.0 which has a soliton centered at max[g(z)]. The solid curve and the points correspond, respectively, 
to the analytical prediction and the numerical result for the imaginary eigenfrequency. 



